{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4ad27792",
   "metadata": {},
   "source": [
    "### test read data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50aec582",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from icecream import ic\n",
    "from pathlib import Path\n",
    "import pathlib\n",
    "\n",
    "_CWD = Path().cwd()\n",
    "_PRJROOT = _CWD.parent.parent\n",
    "_G4OUT = _PRJROOT / \"G4\" / \"build\" / \"Release\"\n",
    "data_path = _G4OUT / \"LDose_ZAll_275.00MeV.npz\"\n",
    "linear_data = np.load(data_path)\n",
    "## Z bin; Phi bin; R Bin; 右边的指标变化最快\n",
    "ix0, ix1, ix2, val, val2, entry = linear_data[\"idx0\"], linear_data[\"idx1\"], linear_data[\"idx2\"], linear_data[\"value\"], linear_data[\"val2\"], linear_data[\"entry\"]\n",
    "##xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\n",
    "Zbin: int = 600\n",
    "Zmax: float = 300\n",
    "Zmin: float = -300\n",
    "## 在 radial 方向, 分成 200 layers; 角度方向不分sector, 只有一整圈\n",
    "Rbin: int = 200\n",
    "Rmax: float = 200\n",
    "Rmin: float = 0\n",
    "##\n",
    "Rmax_fit: float = Rmax\n",
    "Rmin_fit: float = Rmin\n",
    "###=========================================\n",
    "ZDiff: float = Zmax - Zmin\n",
    "Zbwth: float = ZDiff / Zbin  # z bin width\n",
    "Zaxis = np.arange(Zmin + Zbwth / 2.0, Zmax, Zbwth)\n",
    "###=========================================\n",
    "RDiff: float = Rmax - Rmin\n",
    "Rbwth: float = RDiff / Rbin\n",
    "Raxis = np.arange(Rmin + Rbwth / 2.0, Rmax, Rbwth)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "312a40cd",
   "metadata": {},
   "source": [
    "### test data layout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a186769",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = _G4OUT / \"LDose_ZAll_275.00MeV.npz\"\n",
    "iData = np.load(data_path)\n",
    "## Z bin; Phi bin; R Bin; 右边的指标变化最快\n",
    "ix0, ix1, ix2, val, val2, entry = iData[\"idx0\"], iData[\"idx1\"], iData[\"idx2\"], iData[\"value\"], iData[\"val2\"], iData[\"entry\"]\n",
    "\n",
    "## G4 dump 结果安装 fortran 序; # iZ, iPHI, iR\n",
    "seq = 205  ## 对应 csv 209\n",
    "ic(ix0[seq], ix1[seq], ix2[seq], val[seq], val2[seq])\n",
    "\n",
    "## 对应 csv 132802\n",
    "# 331,1,198,1.23936351494257e-13,1.536021922170802e-26,1\n",
    "seq = 132798\n",
    "ic(ix0[seq], ix1[seq], ix2[seq], val[seq], val2[seq])\n",
    "\n",
    "ic(val.shape)\n",
    "iData_rshpd = val.reshape((Zbin, Rbin), order=\"C\")\n",
    "iData_sum = iData_rshpd.sum(axis=1)\n",
    "\n",
    "ic(iData_rshpd.shape)\n",
    "ic(iData_rshpd.flags[\"F_CONTIGUOUS\"])\n",
    "ic(iData_sum.shape)\n",
    "plt.plot(iData_sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64d1b0e0",
   "metadata": {},
   "source": [
    "## test read data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f8f7908",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = _G4OUT / \"IDD_275.00MeV.npz\"\n",
    "iData = np.load(data_path)\n",
    "## Z bin; Phi bin; R Bin; 右边的指标变化最快\n",
    "ix0, ix1, ix2, val, val2, entry = iData[\"idx0\"], iData[\"idx1\"], iData[\"idx2\"], iData[\"value\"], iData[\"val2\"], iData[\"entry\"]\n",
    "\n",
    "## G4 dump 结果安装 fortran 序; # iZ, iPHI, iR\n",
    "iIDD = val\n",
    "\n",
    "ic(iIDD.shape)\n",
    "plt.plot(iIDD)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e27522b",
   "metadata": {},
   "source": [
    "## test functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04d3e8fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "###00000000000000000000000000000000000000000000000\n",
    "ic(val.shape)\n",
    "ic(val.flags[\"C_CONTIGUOUS\"])\n",
    "## G4 dump 结果排布次序为 C order, last fastest; # iZ, iPHI, iR\n",
    "data_2D = val.reshape((Zbin, Rbin), order=\"C\")\n",
    "data_intgrd = data_2D.sum(axis=1)  ## dicrete integrad\n",
    "ic(data_2D.shape)\n",
    "\n",
    "##000000000000000000000000000000000000000000000000 N-高斯拟合设置\n",
    "n_gauss = 5\n",
    "# fit results 存储所有深度的拟合结果\n",
    "fit_rslts = np.zeros((Zbin, 2 * n_gauss + 1))\n",
    "\n",
    "# 剂量最大位置\n",
    "idz0 = np.argmax(data_intgrd)\n",
    "peak_idd = data_intgrd[idz0]\n",
    "ic(idz0)\n",
    "ic(peak_idd)\n",
    "\n",
    "\n",
    "z_depth_index = int(0.85 * idz0)  # 选一个peak之前的位置\n",
    "print(\"Zaxis[{}]: {}\".format(z_depth_index, Zaxis[z_depth_index]))\n",
    "\n",
    "\n",
    "def normalize_cylinder(data1D, binwidth):\n",
    "    sum = np.sum(data1D)\n",
    "    if sum == 0:\n",
    "        return np.zeros_like(data1D)\n",
    "    rs = np.arange(len(data1D)).astype(float)  # 创建一个索引数组\n",
    "    rs *= binwidth\n",
    "    # ic(r)\n",
    "    areas = np.pi * (1 + 2 * rs)\n",
    "    areas[0] = np.pi * np.pow(binwidth, 2.0)\n",
    "    return data1D / areas / sum\n",
    "\n",
    "\n",
    "## 归一化数据\n",
    "profile = data_2D[z_depth_index, :]\n",
    "ic(profile.shape)\n",
    "ic(Rbwth)\n",
    "profile = normalize_cylinder(profile, Rbwth)\n",
    "plt.plot(profile[:25])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "510e7b8f",
   "metadata": {},
   "source": [
    "## test dataclass GaussModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4f252cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyauto.fit.multigauss import cvt_params_to_stctd, cvt_stctd_to_params\n",
    "import numpy as np\n",
    "from icecream import ic\n",
    "# from rich import print\n",
    "## Example usage demonstrating the parameter bridge\n",
    "\n",
    "# Example raw parameter array: [normFactor, w2, w3, s1, s2, s3]\n",
    "raw_params = np.array([1.2, 0.3, 0.2, 1.5, 2.0, 3.0])\n",
    "n_gauss = 3\n",
    "with_background = False\n",
    "\n",
    "# Convert to structured format\n",
    "structured_params = cvt_params_to_stctd(raw_params, n_gauss, with_background)\n",
    "\n",
    "ic(raw_params)\n",
    "ic(structured_params.norm_factor)\n",
    "ic(structured_params.weights)\n",
    "ic(structured_params.sigmas)\n",
    "ic(structured_params.w0)\n",
    "\n",
    "## Convert back to raw format\n",
    "print()\n",
    "raw_converted = cvt_stctd_to_params(structured_params, n_gauss, with_background)\n",
    "## 比较一下\n",
    "print(f\"raw_converted: {raw_converted}\")\n",
    "print(f\"raw_params   : {raw_params}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "078b0156",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pyauto",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
